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PREFACE 


This is a final report summarizing the work performed under Phase HI (third year) 
of a three phase (three year) research study entitled, "Orbital Maneuvering Vehicle 
(OMV) Thrust Chamber Performance". The study has been performed by CFD 
Research Corporation under NASA MSFC Contract: NAS8-37196. Mr. Klaus Gross 
of the EP55 branch. Propulsion Lab, was the project monitor. 

The authors wish to acknowledge the technical guidance of Mr. Klaus Gross of 
MSFC and Dr. Ashok Singhal of CFD Research Corporation. The authors also wish 
to acknowledge the support of TRW, Inc., in providing the VTE engine geometry 
details and some test data. The guidance of Dr. R. Avva of CFD Research 
Corporation in implementing the radiation model is appreciated. The effort of Ms. 
J.L. Swann in preparing this document is also appreciated. 
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PROJECT SUMMARY 


This is the final report documenting the progress made during Phase III of a 3 year 
program to develop a comprehensive modeling methodology for evaluating the 
performance of the Variable Thrust Engine (VTE). Interim reports have been 
delivered at the end of Phases I and II. 

Design and development of liquid rocket engines is currently based on overall 
performance characteristics at design operating conditions. The VTE is expected to 
operate over a range of power levels (from 10% to 100%) without sacrificing much 
efficiency /performance. The design of such engines requires a good understanding 
of the interactions between various complex physical phenomena such as 
atomization, spray dynamics, vaporization, convective/radiative heat transfer, 
turbulent mixing, and chemical reactions. The objective of this study was to 
develop and incorporate models for the above phenomena into an existing CFD 
code to facilitate detailed simulations of the VTE. The major highlights of this study 
are: 

1. The development of a non-orthogonal, body-fitted coordinate (BFC) 
model to reproduce the geometry of the thrust chamber /nozzle, the 
pintle, and the heat shield. 

2. The development and incorporation of a 2-liquid, two-phase, Eulerian- 
Lagrangian spray model into the CFD code. The spray model is fully 
coupled with the gas model. The droplet sizes are estimated from a 
linear surface wave model for impinging jets. 

3. The development of a chemical equilibrium model to simulate the 
hypergolic reaction between MMH and NTO. Thirteen species were 
considered in this model. The equilibrium chemistry model can also 
be used for other propellant combinations such as Hydrogen and 
Oxygen. An instantaneous reaction model and a finite-rate chemistry 
model are also available in the code. 
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4. The development of a discrete-ordinate radiation model to predict the 
effects of radiative heat transfer on the performance of the VTE and 
study of the effects of radiation heat shield. 

5. A series of simulations were performed to estimate the effects of drop 
sizes, spray dynamics, chemistry, radiative heat transfer, on the VTE 
performance. 

6. The simulations were performed at three different power levels 
namely 10%, 50%, and 100%. 

The results indicate that the use of a gas-gas model to analyze the VTE implies an 
inappropriate approximation and the computed flow field is not realistic. It is 
necessary to use a 2-liquid spray model to simulate the droplet dynamics in the 
chamber. The simulations show that smaller droplets result in faster evaporation 
and better mixing, leading to improved engine performance. Larger droplets 
penetrate deep into the flow and impinge on the splash plate giving rise to an 
uneven vaporization distribution. The use of a chemical equilibrium model was 
necessary to bring predictions in better agreement with the data. Simplified finite 
rate models are unable to account for the generation and production of a number of 
species encountered in the VTE chamber. The results also indicate that radiative 
heat transfer is a small fraction of the convective /conductive heat transfer within 
the VTE and at the chamber /nozzle wall. 

This report presents detailed descriptions of the models, simulations, and discusses 
the scope of desirable future work. 
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NOMENCLATURE 


Section 3.1 


Hj chemical potential 

Hj° standard state Gibbs function 

Pj partial pressure of species 'j' 

aij number of atoms in 'i' and 'j' , 

R universal gas constant 

«/ number of moles of species 'j' 

Cj molar concentration of species 'j' 

Section 3.2 


Cl direction of propagation of radiation beam 

I radiation beam intensity 

k absorption coefficient of the gas medium 
cr scattering coefficient of the gas medium 
<t> phase function for gas scattering 

e wall emissivity 

p wall reflectivity 

p, 77 direction cosines of the radiation beam in x, y, r, directions 

( 0 m weighting function for the mth ordinate 
lb black body radiation 

x , y' local grid coordinates 

/ Jacobian of coordinate transformation 

F grid transformation vectors 

8 angular redistribution function 

V volume of the computational cell 

cc weighting function for the numerical scheme 
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Section 3.3 


m d mass of the droplet 

v droplet velocity 

U gas velocity 

Ad droplet surface area 

Vd volume of the droplet 

g gravity vector 

Cd droplet drag coefficient 

Re Reynolds number 

S h Sherwood number 

d diameter of the droplet 

D mass diffusion coefficient 

x v ,x a vapor mass fraction at droplet surface and free stream 
Sc Schmidt number 

Cd specific heat of the droplet 

Td,Tg droplet and gas temperatures 

q sensible heat transferred to droplet 
k thermal conductivity of the gas 
Nu Nusselt number 

Pr Prandtl number 

h v enthalpy of vaporization of the droplet 
y k mass fraction of species 'k' 

Section 3.4 


<y standard deviation of turbulence probability distribution function 
k turbulent kinetic energy 

t e turbulent eddy lifetime 

t r eddy transit time 

( e eddy characteristic length scale 

r droplet relaxation time 
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Section 3.5 


H thickness of the liquid sheet 

0) growth rate of the instability wave 

<7 surface tension 

k wavenumber of the instability wave 

Abbreviations 


MMH 

MonoMethyl Hyradazine 

NTO 

Nitrogen Tetra-oxide 

OMV 

Orbital Maneuvering Vehicle 

VTE 

Variable Thrust Engine 

BFC 

Body Fitted Coordinate 

SMD 

Sauter Mean Diameter 

PISO 

Pressure Implicit Split Operator 

SIMPLE 

Semi-Implicit Method for Pressure-Linked Equations 

MUSCL 

Monotonic Upwind Scheme for Conservation Law 

TDK 

Two-Dimensional Kinetics 
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INTRODUCTION 


1.1 Motivation and Objective 

The Variable Thrust Engine (VTE), developed by TRW, for the Orbit Maneuvering 
Vehicle (OMV) uses a hypergolic propellant combination of Monomethyl 
Hydrazine (MMH) and Nitrogen Tetroxide (NTO) as fuel and oxidizer, respectively. 
The propellants are pressure fed into the combustion chamber through a single 
pintle injection element. The performance of this engine is dependent on the pintle 
geometry and a number of complex physical phenomena and their mutual 
interactions. The most important among these are: 

1. Atomization of the liquid jets into fine droplets; 

2. The motion of these droplets in the gas field; 

3. Vaporization of the droplets; 

4. Turbulent mixing of the fuel and oxidizer; and 

5. Hypergolic reaction between MMH and NTO. 

Each of the above phenomena by itself poses a considerable challenge to the 
technical community. In a reactive flow field of the kind occurring inside the VTE, 
the mutual interactions between these physical processes tend to further complicate 
the analysis. 

The reliable and optimum design of these engines requires a good understanding of 
the internal reactive flow field. Presently, the development of such engines is based 
on overall operation and performance characteristics. Some measurements are 
available for quantities, such as temperature and pressure at a few accessible 
locations. 

The objective of this project is to develop a comprehensive mathematical modeling 
methodology to analyze the flow field within the VTE. This model can be used to 
obtain detailed information on the internal distributions of temperature, pressure, 
flow direction, and state/chemical composition. This information can then be used 
to optimize design parameters and thus improve the performance of the engine. 
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1.2 Project Tasks 


This project is a three-year, three-phase effort to develop the analytical tool. Each 
phase is approximately of one year duration. The major tasks achieved during this 
study are: 

Phase I 

1. Study of OMV/VTE Design and Physical Processes; 

2. Code Adaptation; 

3. Preliminary Computational Analysis; and 

4. Formulation of the Impinging Injector Model. 

A Phase I Interim Report 1 was prepared to document the results of the above 
mentioned activity. 

Phase II 

1. Development of a Spray Model with Two-Liquid, Two Phase capability; 

2. Implementation of Models for Finite Rate Evaporation and Heat Transfer; 

3. Development of a Chemical Equilibrium Model to Simulate the 
Hypergolic Reaction between MMH and NTO; 

4. Analysis of the Flow Field within the Heat Shield Area; 

5. Stability analysis to estimate drop sizes; and 

6. Analysis for Different Power Levels. 

The details of the above mentioned tasks were documented in a Phase II Interim 
Report 2 . The results obtained during these two phases have been presented in a 
JANNAF Combustion Meeting 3 and in an AIAA meeting. 4 

Phase III 

1. Review of existing radiation heat transfer models and selection of a 
model for incorporation into the code; 

2. Development of a turbulent dispersion model for the droplets; 

3. Performing the OMV/VTE thrust chamber calculations using the new 
capabilities for different power levels (100%, 50%, and 10%); 

4. Investigation of the effect of radiation on the flow field; and 
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5. Assessment of current capabilities, identification of limitations, and 
recommendations for future work. 

The details of the above mentioned tasks are described in this report. Section 2 
describes the numerical approach and computational grid generation. The physical 
models are presented in Section 3. The validation of radiation and chemical 
equilibrium models is present in Section 4. The results of the simulations are 
discussed in Section 5. The conclusions from this study are presented in Section 6. 
Section 6 contains an assessment of current capabilities and recommendations for 
future work. 

1.3 Accomplishments of the Project 

The CFD code, REFLEQS, used in the first two Phases of this studyb2 was based on 
the staggered grid approach of the finite-volume method. This traditional staggered 
grid approach is a well known and well established methodology. However, recent 
advances 5 in the colocated grid approach makes it more attractive and more suitable 
for flows in complex geometries. The staggered grid approach involves two control 
volumes: one for the scalar variables and the other for velocities. All the scalar 
variables are stored at the cell center and the velocities are stored at cell faces. In the 
colocated grid approach, all the variables share the same control volume and all the 
variables are stored at cell center. The two approaches have comparable accuracy 
and convergence rates. The colocated grid approach, however, has several 
advantages: 

a. The convection contribution to the coefficients in the discretized 
equations is same for all the variables. This increases computational 
efficiency for higher order schemes and for 3-D geometries. 

b. For complex geometries, cartesian velocity components can be used in 
conjunction with non-orthogonal coordinates, yielding simpler 
equations than when grid-oriented velocity components are employed. 

c. The colocated approach can be extended to multigrid approach. 
Recirculating flows often require very fine grids to achieve accurate 
solutions. The multigrid approach is a way to accomplish this at a 
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reasonable computing cost. Thus, the colocated approach has 
advantages when these techniques are used, 
d. The treatment of boundary conditions are relatively simple for 
colocated approach because of two reasons: no half cells are needed and 
the control volume faces coincide with boundaries. 

The colocated grid approach has been implemented in the REFLEQS code during the 
third Phase of this study. This represents a major enhancement of the code 
capabilities. 

The second major accomplishment is the development of a radiation model for BFC 
geometries. A literature survey was performed to select a suitable model for 
predicting radiative heat transfer in rocket engines. Most of the work reported in 
the literature on radiative heat transfer deal with orthogonal geometries and some 
of the methods reported cannot be extended to complex geometries. The discrete- 
ordinate radiation model, however, was found to be very suitable for extension to 
complex geometries. Also, the discrete-ordinate method can be used for strongly 
participating media such as the combustion gases in the VTE engine. The 
methodology for extending this method to BFC geometries is original and will be 
published in the literature in the near future. With the implementation of BFC 
radiation model in the code, the radiative heat transfer in rocket nozzles can be 
predicted more accurately. In addition to the spray model development, the 
equilibrium chemistry model developed during the second phase of the project has 
been validated with ODE and TDE predictions. 

Finally, several enhancements were made in the spray model. The evaporation 
rates are predicted more accurately by including the effect of liquid saturation 
pressure variation with temperature. The droplet trajectory calculation for BFC 
geometries are made accurate by tracking the droplets in the physical plane. A 
droplet turbulent dispersion model has been included in the spray model. Various 
drop size distribution functions can be used in specifying the droplet initial 
conditions. 

The design and development of VTE engine was mostly based on overall operation 
and performance characteristics. This study has provided a detailed knowledge of 
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what exactly happens in the combustion chamber and nozzle of VTE. This 
information can be used to optimize design parameters and thus improve the 
performance of engine. 
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2 . 


COMPUTATIONAL APPROACH 


This section describes the numerical approach, the geometry and flow conditions 
used in the calculations, and the code adaptation that was done for this study. 

2.1 Numerical Approach 

An advanced Computational Fluid Dynamics (CFD) code, REFLEQS 6 , was used for 
the numerical simulation. The salient features of REFLEQS are given below. 

1. Finite-volume approach of solving Favre-averaged Navier-Stokes 
equations. 

2. Cartesian, polar, and non-orthogonal Body-Fitted Coordinates (BFC). 

3. Colocated (non-staggered) grid. 

4. Strong conservative form of momentum equations with Cartesian 
components as dependent variables. 

5. Pressure-based methodology for incompressible and compressible 
flows. 

6. Fully implicit and conservative formulation. 

7. Pressure-based solution algorithms, including SIMPLE, a variant of 
SIMPLEC, and non-iterative PISO algorithm for all flow speeds. 

8. Four differencing schemes: upwind, central (with damping terms), 

MUSCL, and third-order Osher-Chakravarthy. 

9. Steady-state and transient (second-order time-accurate) analysis 

capability. 

10. Symmetric Whole-Field linear equation solver and Conjugate 

Gradient Squared solver. 

11. Stationary and rotating grid systems. 

12. Multi-component flows with heat and mass transfer, and chemical 
reactions. 

13. Eulerian-Lagrangian approach for two-phase flows. 

14. Range of turbulence models including k-e model, Baldwin-Lomax 
model and the Low-Re model. 

15. Discrete ordinate thermal radiation model. 
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Table 2-1 provides an updated list of features and capabilities of the REFLEQS code. 


2.2 Geometry and Computational Grid 

Figure 2-la presents the overall configuration of a (VTE) in which the main 
combustor and nozzle are shielded by a conical heat shield. The internal surface of 
the heat shield is covered with a reflective insulation blanket for the main chamber 
and nozzle. During the engine operation in vacuum environment, the volume 
between the nozzle wall and the shield surface may be filled with recirculatory hot 
combustion products. 

The VTE thrust chamber, shown schematically in Figure 2-lb consists of a 
cylindrical nozzle and an injector head. A single impinging injector is mounted 
axially at the center of the injector head. The engine is supplied with hypergolic 
propellants MMH (monomethyl-hydrazine CH 3 NHNH 2 ) as fuel and NTO 
(Nitrogen Tetroxide (N 2 O 4 )) as oxidizer. The chamber external wall is cooled by 
radiation while the propellant spray flow produces a fuel rich environment adjacent 
to the inside wall and enhances its cooling. 

The configuration of the OMV/VTE thrust chamber was obtained from detailed 
design drawings. The data from the drawings were used to generate a Body Fitted 
Coordinate (BFC) computational grid shown in Figure 2 - 2 a. Figure 2 - 2 b shows the 
magnified view of the combustor section. Figure 2-2c shows the magnified view of 
the pintle geometry. The computational grid in Figure 2-2a has 100 x 30 cells in the 
longitudinal and radial directions, respectively. Figure 2-3 shows the grid for the 
nozzle with the heat shield. The computational domain has been extended 
downstream of the nozzle exit to resolve the plume. 
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Figure 2-1 (a). Schematic of VTE Combustion Chamber and Nozzle 



Figure 2-1 (b). Schematic of VTE Thrust Chamber 
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(a) Computational Grid for the OMV/VTE Combustor and Nozzle 



(b) Magnified View of the Combustor Grid 



(c) Magnified View of the Pintle Grid 
Figure 2-2. Computational Grid for the OMV/VTE Engine 
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Figure 2-3. Computational Grid for the OMV/VTE Engine with Heat Shield 

2.3 Boundary and Flow Conditions 

A constant mass flow rate is prescribed at the oxidizer and fuel inlets of the VTE 
combustor/nozzle. The density of the incoming propellants is calculated using the 
reference chamber pressure an inlet temperature of 300°K. For liquid injection, the 
drop sizes are prescribed according to a presumed distribution function. The 
droplets are initially located approximately around the point of impact of the two 
liquid jets. The droplets are injected at a temperature of 300 °K. Since the flow at the 
exit is supersonic, the pressure at the exit plane is extrapolated from the cell center 
pressure. The walls of the VTE combustor/nozzle are treated as adiabatic. A 
standard k-e model is used to model the turbulent flow field. A chemical 
equilibrium model with thirteen species is used to simulate the hypergolic reaction 
between MMH and NTO. 

Physical properties for the propellants (MMH and NTO) are obtained from data 
provided by TRW. Some of the properties of these propellants are listed in Table 2- 
2. The mass flow rates of the propellants and the injection areas are obtained for 
operation at the 100% power level. For reference, the conditions used in the 
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computations are summarized below in Table 2-3. The mass flow rates and chamber 
pressure are proportionately reduced for the lower power levels. 


Table 2-2. Physical Properties of MMH and NTO 


Property 

MMH 

NTO 

Density, kg/m 3 

870.0 

1433.0 

Boiling Point, °K 

360 

300 

Saturation Pressure, N/m 2 

6.6 x 10 3 

1 .2 x 10 5 

Heat of Vaporization, KJ/kg 

875.0 

414.0 

Heat Capacity, J/kg°K 

3000.0 

1580.0 

Surface Tension, N/m 

33. 9x 10 * 3 

25.1 X 10' 3 


IMS/iti-i' 


Table 2-3. Flow and Reaction Parameters 


NO. 

Power 

Level 

Flow Rate 
(kg/s) 

MMH NTO 

Pressure 

(psia) 

1 

100% 

.089 

.146 

100 

2 

50% 

.0445 

.073 

50 

3 

10% 

.0089 

.0146 

10 


407 fv 3 12-2 


2.4 Code Modification and Adaptation 

The code was modified during this study to include the following models: 

a. a 2-liquid, two-phase Eulerian-Lagrangian spray model; 

b. a droplet turbulent dispersion model; 

c. a chemical equilibrium model to simulate the hypergolic reaction 
between MMH and NTO; 

d. a thermal radiation heat transfer model; and 

e. a linear stability analysis model to predict drop sizes. 

The details of each of these models are discussed in the following sections. 
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3. 


PHYSICAL MODELS 


This section presents detailed technical descriptions of the physical models 
developed /adopted during this study. 

3.1 Chemical Equilibrium Model 

In the previous phases of this project, the hypergolic reaction between MMH and 
NTO was modeled using single global step instantaneous and finite-rate chemistry 
models. It was shown that these models overpredict the flame temperature. The 
reason for this is that the actual reaction between MMH and NTO involves about 
200 species and several reaction steps. Detailed modeling of all the species and 
reaction steps would be computationally very expensive. It was felt that a chemical 
equilibrium model would be more appropriate in terms of providing reasonable 
accuracy as well as computationally efficient. After reviewing the literature, the 
equilibrium model based on the element potential method was incorporated into 
the REFLEQS code. This model is very general and thus can be applied to any 
propellant combinations provided that important species in the reactions are 
included. This model is very similar to the equilibrium model in the ODE code. 

The chemical equilibrium model proposed by Morgan and MeintjesS is described 
below. This method is based on calculating the element potentials or the element 
activities and is an indirect approach of minimizing the Gibbs function of the 
system. 

The chemical potential of any species can be written as: 


where fij = chemical potential, [i° = standard state Gibbs function, p 0 = atmospheric 
pressure and pj = partial pressure of species . 



(3.1) 



(3.2) 
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where tijln is the mole fraction of species 'f , p is the total pressure of the mixture 
and p 0 is the atmospheric pressure. 


The chemical potential "p" can be written as: 


i=l 


(3.3) 


where //,• is the chemical potential of atomic element i and fl,y is the number of atoms 
of T in species Substituting into Equation (3.1), 


1 


i 


RT ^ 

nt= 1 


a ; 




LL 

RT 


+ /»p-2- 


(3.4) 


From the equation of state. 


Therefore, 


p _RT 
n~~V 


(3.5) 





(3.6) 


where rij/V = molar concentration = Cj. 


Rewriting Equation (3.6), 


In Cj = - 


RT 



(3.7) 
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where Cj = rij/v and A; = fijRT. 


Equation (3.7) can be written as: 


In c- = In Rj + X a { j In Z, 

i= 1 


where 


p P ° 

R rla ex r 



and 


Z, = (an element variable) 


Equation (3.8) is simplified as: 


(3.8) 


(3.9) 


i=l 


(3.10) 


Since atoms of each element are conserved, we can write: 

I 

X a ij c j = b i ,i = l , .../ 

/=* 

where fc, is the total number of atoms of element "i" in all species. 


(3.11) 


Equations (3.55) and (3.56) form the basis of this formulation. Substituting Equation 
(3.55) into (3.56) we get 


£Lr ; n (Z m ) 4 -b,. = 0 , / = !,.../ 

,= I \ m = l ) 


(3.12) 
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Thus we have T non-linear algebraic equations for the T unknown values of Z,. 
The values of cy are calculated from the solved values of Z, (using Equation 3.10). 
The temperature is calculated from cy and then the new values of Ry are calculated. 
The system is solved again for new values of Z; and so on until convergence is 
achieved. 

3.2 Radiation Model 

A number of numerical techniques 8 ^ are available for solving the equation 
governing the transfer of thermal radiation. The most widely used among them are 
listed below. 


a. Hottel's Zone method 10 

b. Flux method 11 

c. P-N approximation method^ 

d. Monte-Carlo method 18 

e. Discrete-transfer method 14 

f. Discrete-ordinate method 18 -^ 

Hottel's zone method was the first general numerical procedure developed to 
handle the radiative transfer equation. It has been widely used to estimate the 
radiation transfer in the absence of detailed knowledge of the flow and reaction. 
This method involves the use of exchange factors, which needs to be worked out in 
advance for complex geometries. The disadvantage of this method is that it is very 
uneconomical and not suited for complex configurations with strongly participating 
media. The flux methods involve a differential approximation to the radiative 
transfer equation. This method is computationally efficient and has, therefore, been 
widely used in general combustion problems. However, this method suffers from 
inadequate accuracy and the difficulty to extend to complex geometries. Compared 
to this method, the P-N approximation method is more accurate. 

Among all the existing numerical treatments for radiative heat transfer, the Monte- 
Carlo method is the most accurate method. In this method, the randomly chosen 
energy releases are tracked through the combustion chamber for their lifetimes. 
Unfortunately, this method is computationally expensive. The discrete-transfer 


16 


4075/3 


method, on the other hand, is economical and can be applied in a straightforward 
manner to complex geometries. This method, however, cannot treat the scattering 
by the medium, containing droplets and soot particles, and hence is not suitable for 
combustor applications. 

After comprehensive review of all the methods, it has been found that the discrete- 
ordinate model is well suited for radiation prediction in rocket combustion 
chambers. There are several advantages of this method. Firstly, the discrete- 
ordinate method can be easily extended to polar and complex geometries. Secondly, 
the evaluation of the in-scattering term is relatively simple. Also, the boundary 
conditions can be imposed more accurately at inlets and exits. Finally, this method 
is accurate and can be easily coupled with CFD codes, based on the finite volume 
approach. The governing equations and the solution procedure for this method are 
described below. 


3.2.1 Governing Equations 

The integro-differential radiative heat transfer equation for an emitting-absorbing 
and scattering gray medium can be written asi3 

(£2 . V) J(r, «}=-{* + o)/(r, £2) + k I b (r) + ^ J[._ „ l{r , £2") 4> (£2' -> £2 ) i £2' m3) 

where 12 is the direction of propagation of the radiation beam, I is the radiation 
intensity, which is a function of both position (r) and direction (12), k and <7 are the 
absorption and scattering coefficients respectively. Ib is the intensity of black body 
radiation at the temperature of the medium, (f> is the phase function of the energy 
transfer from the incoming 12' direction to the outgoing direction 12. The term on 

the left hand side represents the gradient of the intensity in the specified direction 12. 
The three terms on the right hand side represent the changes in intensity due to 
absorption and out-scattering, emission and in-scattering respectively. 


The boundary condition for solving the above equation may be written as 


/(r, Q) = el b (r) + — f \n-Q.\l(r,Q]dQ. 

ft Jn SI <0 


(3.14) 
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where I is the intensity of radiant energy, leaving a surface at a boundary location, e 

is the surface emissivity, p is the surface reflectivity, and n is the unit normal vector 
at the boundary location. 

In the discrete-ordinate method. Equations (3.13) and (3.14) are replaced by a discrete 
set of equations for a finite number of ordinate directions. The integral terms on the 
right hand side of the above equations are approximated by summation over each 
ordinate. The discrete-ordinate equations may then be written as, 

^ W iti §mm An ' ' m ~ ^ (3 ^5) 

m 

For cylindrical polar geometries, the discrete-ordinate equation may be written as. 


BL 


BL 


77 ="(*+<>) An + ^ + ^m 




4 n 


Am ^ ( r An) 2 ^ nJmj m w ; r ° V at 1 \a 


(3.16) 


where the second term on the left hand side represents the angular intensity 
redistribution. It accounts for the change in the direction cosines as a beam travels 
in a straight line. 

The different boundary conditions, required to solve this equation, are given as: 

Wall Boundary 


/ m = e/ b + 7 2 w m \ for west boundaries 

vL<o 


I — £ /l + — 
m b n 


^ / / ' 

2j w m | x m I m ;for east boundaries 
y- m >o 
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C<o 



C>o 


(3.17) 


Symmetry Boundary 


l m = l m ; \i m = - \i m for east and west boundaries; 

I m = / m ' ; £, m =-^, m for south and north boundaries 


(3.18) 


Inlet and Exit Boundary 


(3.19) 


In the above equations, m and m' denote the outgoing and incoming directions, 

respectively. For a particular direction, denoted by m, the values /i and £ are the 
direction cosines along the x and y directions. These equations represent m coupled 
partial differential equations for the m intensities, l m . In the present code, only the 
S4 approximation, i.e. 12 ordinate directions, is considered. Table 3-1 lists the values 
of the direction cosines and the associated weights. 

Equation (3.15) can be cast in a fully conservative form for a general BFC coordinate 
system as 


[n m L / F„ + /_ / F,„l +2, rn_ L I F„ + F 2y ] 



(3.20) 


where x and y' are the local grid coordinate directions, / is the Jacobian of the 
coordinate transformation and 
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( 3 . 21 ) 


Table 3-1. Direction Cosines and the Associated Weights for the S 4 Approximation 



4075/3 t3- 1 


Equation (3.16) can be integrated over the control volume, shown in Figure 3-1, to 
obtain 


+ ^m n ^m s ) + [^n ^s) 


[Ym + 1/2 l P, m+ 1/2 “ Ym - 1/2 ^P,w - 1/ 2] 




= -(fc + a)V/ mP + fcV/^ + V^I«; m '< ^ m I, 


mp 


(3.22) 
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where the coefficients y are evaluated from the recurrence relation 


^m + l/2 m-1/2 (3.23) 

In the above equation A denotes the cell face area and V denotes the cell volume 
and subscripts w, e, s, n, P denote the west, east, south, north, and cell center, 
respectively. Figure 3-1 shows a sketch of the control volume. 



Figure 3-1. Sketch of the Control Volume 

The in-scattering term on the right hand side of this equation contains the phase 

function 0, which can be prescribed for both isotropic and anisotropic scattering by 
the medium. For a linear anisotropic scattering, the phase function may be written 
as, 


<t>mm = 1 • 0 + a 0 [\i m [i m +ri m r| m ] 


(3.24) 
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where a 0 is the asymmetry factor that lies between -1 and 1. The values -1,0,1 denote 
backward, isotropic and forward scattering, respectively. 

3.2.2 Solution Procedure 

Equations (3.22) involves 5 intensities, out of which 2 intensities are known from 
the boundary conditions. The other cell face intensities can be eliminated in terms 
of the cell center intensity by a spatially weighted approximation. For example, if 
the east and north cell face intensities are unknown, they can be eliminated in terms 
of the cell center intensity by using 


2 m n ~ qJ (3.25a) 

l m t = ^ (3.25b) 

where a is the weighting factor, a-1 yields a upwind differencing scheme and a=2/2 
yields a central differencing or the "diamond differencing" scheme. Substituting 
Equation (3.25) into the integrated discrete-ordinate equation results in a algebraic 
equation for the cell center intensity. 


^m AI m w + ^m BI m s + 5 l- S 2 

j W S 

m p " IV^ + ^As + ccffc + aJV (3.26) 


where 
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A = A e (1 - a) + A w a 
B = A n {l-a) + A s a 
Sj=a kVIfy 

G 

S 2 — oc V w m ty mm I m p 


Equation (3.26) is appropriate, when both the direction cosines are positive and the 
direction of integration proceeds in a direction of increasing x and y. For a negative 
direction cosine, the direction of integration is reversed, and the integration sweep 
is started from the appropriate corner of the domain. Thus, for each set of direction 
cosines, the intensities at all the cell centers are obtained by marching in the 
appropriate direction. This procedure is repeated for all the ordinate directions. The 
in-scattering term is evaluated explicitly by using the previous iteration values, thus 
decoupling the m ordinate equations. 

Assuming local radiative equilibrium, the source term for each computational cell 
i.e. the net radiative heat flux, is given by. 


S = Vk'Ll m io m -V4nkI b 

m m b (3.27) 

m 

This source term is added to the gas phase energy equation. The gas phase equations 
and the radiative transfer equations are solved in an iterative manner. 

3.3 Spray Model 

There are two approaches for modeling sprays: (a) Eulerian-Eulerian; and (b) 
Eulerian-Lagrangian. In the Eulerian-Eulerian approach, the droplet phase is treated 
as a continuum. Thus, the liquid phase equations are similar to the gas phase 
equations. These two sets of equations are solved in a strongly coupled manner. In 
the Eulerian-Lagrangian approach, individual drop parcels are tracked until their 
lifetime in a Lagrangian frame of reference. The gas phase equations are coupled 
with the liquid phase equations through droplet source/sink terms. The major 
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assumption used in this approach is that the volume occupied by the liquid phase is 
negligibly small when compared to the cell volume. This assumption does not hold 
in the dense spray region. Thus, the Eulerian-Eulerian approach is more accurate is 
in this respect. However, the Eulerian-Eulerian approach suffers from other 
difficulties such as boundary condition treatment, smearing of the liquid phase, 
effect of gas phase pressure gradient, etc. Moreover, in typical reacting flow 
problems, the spray region is small when compared to the flow domain. Hence, the 
Eulerian-Eulerian approach is computationally inefficient. Due to these reasons, the 
Eulerian-Lagrangian approach is widely used and it is implemented in the REFLEQS 
code. 

The droplet equations of mass, momentum and energy are solved in a Lagrangian 
frame of reference moving with the droplets. These solutions are used to calculate 
the source/ sink terms for the corresponding gas phase equations. 

The equation of motion for the droplet is written as^ 9 : 


civ V Q v r v Q 

m d-ji = C D P{U-v)\U-v\ — -Vp-V d + m d — e r -m d — r e Q + m d g ( 3 . 28 ) 

where is the mass of the droplet and v = v x e x + v r e r + veee its velocity vector. 
v x , v r and vq are the axial, radial and swirl components of the velocity. Co is the 

drag coefficient and p, U, and p are the density, velocity, and pressure of the gas, 
respectively. Ad is the droplet surface area, and Vd is the droplet volume. For a 
spherical droplet. Ad = tuH and Vd = (7cd 3 /6), where d is the droplet diameter, g is the 
gravity vector, and r is the radial position of the droplet. Equation (3.28) accounts for 
the acceleration/ deceleration of the droplet, due to combined effects of drag with gas, 
local pressure gradients, centrifugal and Coriolis forces, and body forces, such as 
gravity. The centrifugal force term appears only in the radial velocity equation, and 
the coriolis force term appears only in the swirl velocity equation. 

The drag coefficient, Co, for the droplet is based on the local Reynolds number, 
evaluated as: 


24 


4075/3 


(3.29) 


Re - el£_- 

where p is the dynamic viscosity of the gas. The following correlations have been 
found to be valid for different ranges of Reynolds numbers: 


c D = |i 

Re 

for Re < 1 

C D = ^-[ 1+ 0.1 5Re 0 - 687 ] 
Re 

for 1 < Re < 10' 

C D = 0.44 

for Re > 10 3 


Substituting for Cp into Equation (3.28), 


dv 

dt 


p/ 


(U-v). 


vp v i 
— + — 
P d r 


v r v Q 

“7 _e e + s 


(3.30) 


(3.31) 


where pd is the droplet liquid density and /= CoRe/24. The droplet locations are 
determined from the velocities as: 


dR = v 
dt 


(3.32) 


where R is the position vector of the droplet. 


The mass conservation equation for the droplet is given by: 


^^ = -Sh(pD)7td(x v -xJ (3.33) 

dt 

where Sh is the Sherwood number, and D is the mass diffusion coefficient for the 
gas. x v and Xo° are the vapor mass fractions at the droplet surface and the free 
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stream, respectively. In this model, x v is calculated from expressions for the vapor 
pressure as a function of droplet temperature. The Sherwood number is obtained 
from the following expression: 

Sh = 2 + 0.6Re°- 5 S c 0333 (334) 

where S c is the Schmidt number. Since md = pd7a/ 3 /6, the mass equation for the 
droplet is rewritten in terms of its diameter: 


did) 

dt 


-2Sh (pD) 



The energy equation for the droplet is written as: 


(3.35) 


mdCd = i + (3.36) 

where q is the sensible heat transferred to the droplet and Tj is the droplet 
temperature. L is the latent heat of vaporization for the droplet fluid and Q is the 
specific heat of the droplet, q is calculated as: 

q = NuKkd {T g -T d ) (3.3 7) 

where k and T g are the thermal conductivity and temperature of the gas, 
respectively. The Nusselt number, Nu, is obtained from the following correlation: 


Nu = 2 + 0.6Re°- 5 Pr 0 - 333 


(3.38) 


where Pr is the Prandtl number. Substituting the above expressions into Equation 
(3.36), the energy equation is rewritten as: 


dT d _ Tg~ Td Q l 

dt q e 


(3.39) 
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where 


Ql = 


L Sh {pD){x v -x„) 
NuJt 


p^ 2 c d 

e =^Tu T 


(3.40) 


(3.41) 


If the droplet temperature is below its boiling point. Equation (3.35) is used as the 
mass conservation equation. Once the droplet attains its boiling point, the mass 
conservation equation is obtained by setting the left hand side of Equation (3.36) to 
zero, i.e., the droplet temperature cannot exceed its boiling point. Therefore, 


dm_d _ ? 
dt ’ L 


(3.42) 


This equation implies that all the heat transferred from the gas is used to vaporize 
the droplet. 


3.3.1 Calculation of Droplet Source/Sink Terms 

The solutions to the droplet equations provide sources/sinks to the gas phase 
conservation equations. The sources/ sinks are assigned to the particular cell (in the 
Eulerian gas phase) in which the droplet is located. The calculation of the 
source/ sinks is as follows. 


The mass of the droplet is continuously monitored along its trajectory. The 
source/sink term for a cell is given as: 


Am = p d K n X 


(4 -(4 

v win V wout 


(3.43) 
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the summation takes place over all the droplets located in the given cell. The 

subscripts "in" and "out" refer to cell inlet and outlet conditions, and n is the 
droplet number density. Similarly the momentum source/sink term is evaluated 
as: 


AM = p d K n X 


ML -ML 


(3.44) 


This includes the effects of frictional interactions with the gas. The momentum 
source term is a vector. The energy source/ sink term is given as: 


A £ = Amh v - n £ c\ { 

i 


(3.45) 


where h v is the enthalpy of the vapor at the vaporization temperature and q. is the 
heat transfer from the gas to the droplet (calculated according to Equation (3.37)). 

3.3.2 Governing Equations for the Gas Field 

The mass conservation for the flow can be represented by the continuity equation. 
This formulation accounts for a variable density flow field: 


V-pu = S p 


where 


S 


p 


Am 


(3.46) 


(3.47) 


Am is obtained from Equation (3.43), and V is the volume of the cell, in which the 

mass addition takes place, p is the density of the fluid, and u is the velocity vector. 
The mass conservation for individual species may be written as: 
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VpuY k = V{D t VY k )+W k + S t 


(3.48) 


where V and is the mass fraction of specie, k, and W k is a source term. 


Am k 


consisting of reaction rate expressions. A m k is the mass of species k added, due to 
vaporization. 


The momentum conservation equation for the flow is written as: 


V ■ p uu = - Vp - V ■ x + pg + S m 


(3.49) 


S - 

where m V , and t is the stress tensor equal to -p(Vu + Vu T ). A is obtained 


AM 


from Equation (3.44). It is assumed that the gas behaves as a Newtonian fluid. The 
energy conservation equation for the system is given by: 


enthalpy and the thermal conductivity of the gas, respectively. 

3.4 Droplet Dispersion Model 

The dispersion of the droplet, due to turbulence, is modeled using the method of 
Gosman and Ioannides. 18 In this method, an instantaneous gas velocity is used in 
the droplet equations of motion. The instantaneous gas velocity is evaluated from 
the time-averaged gas velocity, the turbulent kinetic energy, and dissipation rates. 
For this purpose, the turbulence is assumed to be isotropic and to possess a Gaussian 



(3.50) 
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probability distribution function in the fluctuating velocity, whose standard 
derivation is given by 


a = 



(3.51) 


where k is the turbulence kinetic energy in the computational cell, where the 
droplet is located. Random sampling of this distribution at appropriate times 
intervals along the droplet trajectory yields the fluctuating velocity of the gas and 
hence the instantaneous gas velocity. 

The sampled fluctuating velocity is applied over a time interval called the 
interaction time, in the droplet trajectory calculation. The interaction time can be 
thought of as the time scale, associated with a turbulent eddy. There are two 
possible events: (a) the droplet moves sufficiently slowly relative to the gas to 

remain within the eddy during its lifetime (f e ); (b) the relative velocity between the 
droplet and gas is sufficient to allow it to traverse the eddy in a transit time (f r ). The 
interaction time scale will be the minimum of the time scales in the above events, 
i.e., 


t ini = min (t e/ t r ) (3.52) 

The eddy and transit time scales can be estimated as follows. 

The eddy characteristic length scale is assumed to be the dissipation length scale. 



(3.53) 


where e is the turbulence dissipation rate. 
The eddy lifetime is then estimated as: 
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(3.54) 



The transit time scale, t r/ is estimated from the following solution of a simplified 
and linearized form of the equations of motion of the droplet 


t r = -i In 



iTlU-vl 


where x is the droplet relaxation time given by 


4 d 

T '3 P ‘'(pC D |U-v|) 


(3.55) 


(3.56) 


In the case when [x|U- v|j > ^ / the above equation has no solution. This implies 

that the droplet is captured by the eddy, in which case tint = t e . 


3.5 Linear Stability Model for Imping in g Tet Atomization 


A preliminary stability analysis was conducted to determine the sizes of the droplets 
that would result from the impact of the liquid jets. Daniel Wiehs 21 has obtained 
the following expression for the stability of a radially moving liquid sheet in 
stagnant air: 


w = 


vk 2 {p L H) 1/2 

2 


-1 + 

, ip c ku 2 - ok 2 )V 12 


v*k 4 p L H 1 


where 


(3.58) 


w = growth rate of the surface wave 
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v 

kinetic viscosity of liquid 

P L 

liquid density 

P G = 

gas density 

H 

liquid sheet thickness 

U 

sheet velocity 

a = 

surface tension 

k 

wave number (m* 1 ) 


The liquid density is estimated as an average of the densities of MMH and NTO. 
The liquid velocity is calculated from the mass flow rate of MMH and the cross- 
sectional inlet area. The initial thickness of the sheet (H) is taken to be twice the size 
of the MMH inlet. This expression is valid only for the low speed breakup of the 
sheet, i.e., drop sizes are of the order of the sheet thickness. Figure 3-2 shows the 
variation of growth rate with wavenumber. 



Figure 3-2. Plot of Growth Rate Versus Wave Number 
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However, the above analysis assumes that the gas or the ambient is stationary. This 
is not the case inside the VTE chamber, where strong recirculation zones (resulting 
from the geometry of the chamber) can result in droplets being stripped away from 
the surface of the liquid sheet. The drop sizes resulting from the high speed breakup 
are considerably smaller. The low speed analysis for this problem predicts a 
maximum drop size of about 300 microns. A mean drop size of 50 microns is used 
in the calculations in order to approximately account for high speed atomization 
phenomena. A lognormal distribution (about the mean size) is used to calculate 
other drop sizes. 
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4. 


MODEL VALIDATION 


4.1 Radiation Model 

The discrete-ordinate radiation model, implemented in the present code, has been 
validated with several benchmark cases, for which exact analytical solutions are 
available. First, the radiative heat transfer in a square enclosure with an absorbing 
medium is considered. The walls of the enclosure are maintained at 0°K and the gas 
in the enclosure is maintained at 500°K. Results are obtained for various optical 
thickness. A 10 X 10 grid and the S4 approximation, i.e. 12 ordinate sets, are used for 
all these cases. Figures 4-la through 4-1 c show the non-dimensional radiative heat 
transfer at the cold walls for three optical thicknesses of 0.1, 1.0, and 10.0. The circles 
represent the numerical solution and the lines show the exact analytical solution 
reported in Lockwood and Shah (1981). As can be seen, the agreement is very good. 

Next, the radiative heat transfer in a square enclosure with a strongly participating 
medium is considered. In this problem, the east, west and north walls are 
maintained at 0°K and the south wall is maintained at 1000°K. All the walls are 
black and the gas scattering coefficient is zero. Figures 4-2a through 4-2c show the 
radiation heat transfer at the south wall for absorption cross-sections of 2.0, 1.0, and 
0.1 respectively. A 21 X 21 grid and the S4 discrete-ordinate approximation are used 
for these calculations. The open symbols show the exact solution reported by 
Razzaque et al. (1983), and the lines show present predictions. Again, the model 
results show good agreement with the exact solution. 

Finally, radiative heat transfer for pure scattering in a square grey enclosure is 
considered. This case is quite important, because the boundary conditions are not 
only a function of the surface emissive power but also of the incident radiant 
energy. The boundary conditions for this case are the same as in the previous case. 
A 25 X 25 uniform grid system is used for this case. Figures 4-3a through 4-3c show 
the surface heat transfer at the hot wall for wall emissivities of 1.0, 0.5, and 0.1 
respectively. The comparison with exact solutions are quite good, except for the low 
surface emissivity case. For this case, the error in the intensity calculation affects the 
boundary condition to a large extent. The results can be improved by refining the 
grid and increasing the number of ordinate sets considered. 
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(a) 



(b) 



x/L 

(C) 

Figure 4-1. Surface Heat Transfer Rates for a Square Enclosure with Cold Black 
Walls and an Absorbing Medium for Various Optical Thicknesses: (a) 
0.1, (b) 1.0 and (c) 10.0. o, - numerical simulation; — , exact solution 
from Lockwood and Shah. 14 
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Figure 4-2. Radiative Heat Transfer at the Hot Wall of the Square Enclosure 
Problem for Various Absorption Cross-Sections: (a) 2.0; (b) 1.0; and (c) 

0.1. o, - exact solution 15 ; — , present prediction 
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Figure 4-3. Radiative Heat Transfer at the Hot Wall of the Square Enclosure 
Problem with Grey Walls and Pure Scattering Medium for Various 
Wall Emissivities: (a) 1.0; (b) 0.5; and (c) 0.1. o, - exact solution 15 ; 

— , present prediction 
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4.2 Chemical Equilibrium Model 


The objective of this analysis was to use the chemical equilibrium model to predict 
species distributions in the SSME combustor /nozzle and to compare the predictions 

with the TDK 22 program (TDE option). Figure 4-4 shows the grid (385 x 50) for the 
SSME combustor /nozzle. The boundary conditions for the problem are: 

a. Inlet: 

u = 320.5 m/s 
p = 19.8 MPa 
T = 3585°K 

Composition: 91.1% H 2 0, 4.4% OH, 3.5% 0.56% 0 2 , 0.25% O, 0.19% H 

Figures 4-5a and 4-5b show the temperature distribution predicted by the TDE and 
the REFLEQS codes. It is observed that the REFLEQS predictions match well with 
the TDE solution. Figures 4-6a, 4-6b and 4-6c show the distribution of species (H 2/ O 2 , 
and H 2 0, respectively) along the nozzle centerline. The symbols are the results of 
the TDE program. It is observed that the chemical equilibrium model in REFLEQS 
predicts correctly the variation in concentration of the above species in the nozzle. 



Figure 4-4. Computational Grid for the SSME Nozzle (385x50) 
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Figure 4-5. Temperature Distribution in the Nozzle 
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Figure 4-6. Species Distribution in the Nozzle 
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5. 


COMPUTATIONAL RESULTS 


In this section, the computational results of VTE performance obtained using the 
REFLEQS code are presented. A number of simulations were performed to assess 
the effects of various physical phenomena occurring in the thrust chamber of the 
VTE. These effects are discussed separately in several sections. 

The first section (5.1) discusses the performance of VTE obtained using the gas-gas 
combustion model. The VTE performance is analyzed for three power levels, 
namely 100%, 50%, 10%. An equilibrium chemistry model with 13 species is used 
for all the cases. 

The second section (5.2) considers the effect of atomization and spray. Calculations 
for the three power levels are obtained using a mean drop size of 50 microns with a 
lognormal distribution of sizes ranging from 25-100 microns. A parametric study on 
the effect of drop sizes is also included in this section. An equilibrium chemistry 
model is used for these cases. The spray model used in obtaining these results is 
described in Section 3.3 

The third section (5.3) presents the effect of radiation on the performance and wall 
cooling requirements of VTE. The S4 approximation (12 ordinate directions) of the 
discrete ordinate radiation model is used for obtaining these results. The 
temperature distributions in VTE chamber /nozzle with and without invoking the 
radiation model are compared. 

Finally, in Section 5.4, the effect of the radiation heat shield on heat transfer 
characteristics is discussed. The heat shield is designed to protect the structure 
surrounding the engine from the radiation. For all these calculations, a solution is 
regard to be converged when the residuals of all the variables have dropped by 5 
orders of magnitude. 

5.1 Gas-Gas Combustion Model 

In the Gas-Gas combustion model, it is assumed that atomization and vaporization 
occurs instantaneously at the injection location. In reality atomization and 
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vaporization occur during a finite time scale. But the Gas-Gas model will be useful 
in providing an understanding of the flow and combustion processes occurring 
inside the VTE thrust chamber. For all the VTE performance calculations, an 
equilibrium chemistry model was used to simulate combustion. In the previous 
phase of this study, instantaneous and finite rate chemistry models with five species 
(MMH, NTO, CO 2 , H 2 O, and N 2 ) were considered for the reaction between MMH 
and NTO. The ODE analysis for this system, however, shows that many more 
species are involved in this reaction which needs to be considered for accurate 
prediction of flame temperature. Therefore, a chemical equilibrium model with 13 
species was used for the present calculations. A detailed formulation of this model 
is described in Section 3.1. 

5.1.1 100% Power Level 

For obtaining the following results, the MMH and NTO inlet temperatures are 
taken to be 300°K. The density of gaseous MMH and NTO are calculated from the 
specified chamber pressure, and their velocities are calculated from the specified 
flow rates. The chamber pressure for 100% power level is specified to be 6.89E5 
N/m 2 (100 psi). For most of the calculations performed, a zero gradient boundary 
condition for temperature is imposed at the walls. For problems involving 
radiation, however, a wall temperature distribution is specified. Figure 5-1 shows 
the Mach number contours in the VTE combustion chamber and nozzle. It is 
observed that the maximum Mach number at the exit is about 5.0. It is also be 
observed that the thickness of the boundary layer on the nozzle wall increases with 
downstream distance. Figures 5-2 and 5-3 show the distribution of temperature and 
passive scalar concentration of NTO. The maximum temperature in the chamber is 
seen to be 3300°K which is in good agreement with the calculations of adiabatic 
flame temperature for this mixture. These results are very similar to the results 

reported in the previous report 2 indicating that the accuracy of staggered and 
colocated grid approach are comparable. It should be pointed out, however, that 
calculations performed during the previous phases of this study 1,2 using 
instantaneous and finite-rate chemistry models predicted a much higher 
temperature. 
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Figure 5-1. Mach Number Distribution in VTE Chamber/Nozzle (Gas-Gas Model, 
100% Power Level) 



Figure 5-2. Temperature Distribution in VTE Chamber/Nozzle (Gas-Gas Model, 
100% Power Level) 
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Figure 5-3. Passive Scalar Concentration of NTO (Gas-Gas Model, 100% Power 
Level) 

Figure 5-3 indicates that the reaction between MMH and NTO is complete at the 
location midway between the pintle and the throat. Figure 5-4 shows the velocity 
vectors inside the combustion chamber. Two recirculation zones are observed, a 
smaller one in the recessed portion of the chamber dome, and a larger one near the 
center of the chamber. The presence of this large recirculation zone near the center 
improves the mixing characteristics in the combustion chamber. As a result, the 
temperature downstream of the pintle until the throat is nearly uniform. The 
predicted ISP for this case is about 330 which is about 5% higher than TRW's 
experimental data. This is to be expected because the gas-gas model does not include 
the energy required for evaporating the droplets. 
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VELOCITY PLOTS 
VELMIN 0 000E+00 
VELMAX 1 388E+03 



Figure 5-4. Velocity Vectors in VTE Combustion Chamber (Gas-Gas Model, 100% 
Power Level) 

5.1.2 50% Power Level 

Next, the performance of VTE at 50% power level is considered. The chamber 
pressure for this power level is specified as 3.495E5 N/m 2 (50 psi) and the MMH and 
NTO flow rates are 50% of the full power level. The inlet area for the incoming 
gases is assumed to remain constant for all the power levels. Figures 5-5 and 5-6 
show the Mach number and temperature contours for this case. The maximum 
Mach number at the exit is very close to the maximum Mach number for the 100% 
power level. The maximum temperature in the chamber, however, is about 5% less 
than the maximum temperature for the 100% power level. Except for small 
differences in the magnitudes of Mach number and temperature, the contours for 
both power level cases are very similar. Figure 5-7 shows the velocity vectors in the 
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combustion chamber. Again, the velocity vectors for both power level cases are very 
similar except that the recirculation zones are weaker in this case due to reduced 
mass flow rates. The ISP predicted for this case is slightly smaller than that for the 
previous case mainly due to the off-design conditions. 
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Figure 5-5. Mach Number Distribution (Gas-Gas Model, 50% Power Level) 
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Figure 5-6. Temperature Distribution (Gas-Gas Model, 50% Power Level) 
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VELOCITY PLOTS 
VELMIN 0 000E+00 
VELMAX 1 372E+03 



Figure 5-7. Velocity Vector Plot (Gas-Gas Model, 50% Power Level) 

5.1.3 10% Power Level 

The Mach number and temperature contours for the 10% power level are shown in 
Figures 5-8 and 5-9. The chamber pressure for this case is 10 psi. The propellant 
flow rates are reduced to l/10th of the flow rate at full power level. For the 10% 
power level, the maximum Mach number at the exit is about 4.6 which is 8% less 
than that for the full power level case. The maximum temperature in the 
combustion chamber is 3000°K which is 200°K less than the maximum temperature 
at full power level. Figure 5-10 shows the velocity vectors in the combustion 
chamber. It can be observed that the recirculation zone above the pintle has almost 
disappeared because of the lower velocity of the gases and the recirculation zone at 
the center is very weak. The ISP predicted at this power level is 325. This is slightly 
less when compared to the ISPs at other power levels mainly due to reduced 
efficiency at this condition. The following section presents the results obtained 
using the spray model. 
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Figure 5-8. Mach Number Distribution (Gas-Gas Model, 10% Power Level) 



Figure 5-9. Temperature Distribution (Gas-Gas Model, 10% Power Level) 
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Figure 5-10. Velocity Vector Plot (Gas-Gas Model, 10% Power Level) 

5.2 Liquid-Gas (Spray) Combustion Model 

In this section, the performance of VTE predicted by the spray model is discussed. 
The spray model incorporated in the REFLEQS code is based on the Eulerian- 
Lagrangian methodology of two-phase flows discussed in Section 3.3. For accurate 
prediction of performance, a good estimate of the droplet initial conditions are 
required. The linear stability analysis of impinging jets, discussed in Section 3.5, 
predict a maximum drop size of 300 microns for MMH and NTO liquid jets. This 
analysis is strictly valid for low speed atomization and jets exhausting into still 
ambience. However, in the VTE combustion chamber, the liquid fan formed by the 
impingement of MMH and NTO is surrounded by recirculating and chemically 
reacting gaseous mixture. The breakup process in this situation is expected to be due 
to the high speed atomization mechanism for which the droplet size would very 
small. Also, the secondary atomization process would make the droplet size still 
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smaller. Due to these reasons, the droplet sizes produced by the VTE pintle injector 
are expected to be significantly less than that predicted by the linear stability analysis. 

For the following simulations, a log-normal drop size distribution is specified with a 
mean droplet size (SMD) of 50 microns. The droplet sizes are selected in the 25 to 
100 microns range to fit a log-normal curve with a standard deviation of 0.2. Figure 
5-11 shows the liquid mass fraction as a function of drop size. A parametric study is 
performed to understand sensitivity of predictions to initial drop size distribution. 

5.2.1 100% Power Level 

Figure 5-12 shows the trajectories of droplets for the 100% power level conditions. 
All the drops are released at the point of impact of the two jets. The initial droplet 
velocity is expected to be an average of liquid sheet velocity and gas velocity near the 
interface. In the previous report 2 , a parametric study on initial droplet velocities 
was performed. That study showed that the mixing in the combustion chamber 
depends on the magnitude of the droplet velocity. Higher droplet velocities lead to 
smaller residence times and sufficient mixing is not achieved. On the other hand, 
smaller droplet velocities gives rise to a well mixed zone in the central region of the 
combustor. Based on the results of that study, a mean droplet velocity of 35 m/sec is 
used for all the power level cases with a 5% random variation in velocity is specified 
for each droplet size group. The trajectory plot indicates that the penetration of 
smaller droplets is significantly less than the penetration of larger droplets. It is 
clear that larger droplets impinge on the combustion wall and slide along the wall 
before they evaporate. It is also clear that the spray is very dense near the injection 
point. Figures 5-13 and 5-14 show the Mach number and temperature contours for 
this case. The maximum Mach number at the exit is 4.8 which is slightly less than 
the prediction of gas-gas model. The temperature contours show that there are large 
temperature gradients in the combustion chamber due to the presence of relatively 
cold droplets. In the predictions of gas-gas model, the temperature in the chamber 
was nearly uniform. The hottest zone is in the center of the combustion chamber. 
As a consequence, the head end of the pintle is expected to be very hot. This has 
been proven in the VTE hot fire tests. Later designs of VTE use an ablative material 
cast on the front end of the pintle injector. 
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Figure 5-11. Lognormal Drop Size Distribution for 50p SMD 
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Figure 5-12. Droplet Trajectory Plot (Spray Model, 100% Power Level) 
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Figure 5-13. Mach Number Distribution in VTE Chamber/Nozzle (Spray Model, 
100% Power Level) 



Figure 5-14. Temperature Distribution (Spray Model, 100% Power Level) 
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Figure 5-15 shows the evaporated mass (kg /sec) of NTO in the combustion chamber. 
It is clear than even though some of the droplets impinge on the splash plate, most 
of the evaporation takes place very close to the pintle injector. Figure 5-16 shows 
the velocity vectors in the chamber for this case. As can be seen, there is no 
recirculation zone in the swept back region of the combustor dome. The 
recirculation zone downstream of the pintle is also absent, contrary to the gas-gas 
model prediction. This is expected because the injection velocity of gas is much 
higher than the injection velocity of liquid for the same flow rates. The ISP 
predicted by the spray model is 324 which is slightly less than the ISP predicted by 
the gas-gas model. 



Figure 5-15. Distribution of Evaporated Mass of NTO (Spray Model, 100% Power 
Level) 
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Figure 5-16. Velocity Vectors in VTE Combustion Chamber (Spray Model, 100% 
Power Level) 

5.2.2 50% Power Level 

Figure 5-17 shows the trajectories of droplets for the 50% power level case. The 
spreading angle of the spray is slightly more in this case because of the lower 
chamber pressure. Again the larger droplets do not evaporate before impinging on 
the wall. Figures 5-18 and 5-19 show the Mach number and temperature contours. 
Similar to the predictions of the gas-gas model, the maximum Mach number at the 
exit is slightly lower than that for the 100% power level. The temperature 
distribution for this case is very similar to the distribution for 100% power level. 
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Figure 5-17. Droplet Trajectory Plot (Spray Model, 50% Power Level) 
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Figure 5-18. Mach Number Distribution (Spray Model, 50% Power Level) 
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Figure 5-19. Temperature Distribution (Spray Model, 50% Power Level) 


5.2.3 10% Power Level 


The Mach number and temperature distributions predicted by the spray model for 
10% power level are very similar to the corresponding distributions at other power 
levels. This is very clear from the Mach number distribution shown in Figure 5-20. 


As expected, the ISP predicted for this power level is slightly less than that for other 
power levels. The next section discusses the effect of droplet sizes on the 


performance of VTE. 
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Figure 5-20. Mach Number Distribution (Spray Model, 10% Power Level) 
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5.2.4 Parametric Variation of Droplet Sizes 

The spray model is used to conduct a parametric study to assess the effect of droplet 
size. Two simulations for mean drop sizes of 25 and 100 microns are performed. 
Figure 5-21 shows the drop size distributions used in these two cases. Figure 5-22 
shows the droplet trajectories for the 25 micron size case. It is observed that droplets 
evaporate very quickly and very few droplets impinge on the splash plate. The 
spreading angle of the spray is small compared to the angle for 50 microns size case. 
Figures 5-23 shows the distribution of evaporated NTO mass. It clearly indicates that 
evaporation takes place very close to the pintle. Because evaporation takes place 
quickly, there is sufficient time for mixing and chemical reaction. As a result, the 
temperature distribution in the combustion chamber is more uniform. Figure 5-24 
shows the droplet trajectories for the 100 micron size. As can be seen, more droplets 
impinge on the combustion wall. The large droplets slide along the combustor wall 
and evaporate before reaching the throat. Figure 5-25 shows the contours of 
evaporated NTO mass. Because of the slow evaporation of large drops, combustion 
takes place until the gases reach the throat region. From these simulations, it is 
clear that optimum performance can be achieved when the droplet sizes are in the 
25-50 microns range. 
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5-21. Lognormal Size Distributions for 25g SMD (o) and 100(1 SMD (□) 



Figure 5-22. Droplet Trajectories for 25(1 SMD (100% Power Level) 
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Figure 5-23. NTO Evaporated Mass Distribution for 25 Micron Drop Size (100% 
Power Level) 



Figure 5-24. Drop Trajectories for 100|i SMD (100% Power Level) 
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Figure 5-25. NTO Evaporated Mass Distribution for 100 Micron Drop Size (100% 
Power Level) 


5.3 Effects of Radiation 

In this section, the effect of radiation on flow and heat transfer characteristics in VTE 
is discussed. Calculations were performed by invoking the radiation model, using 
both gas-gas and spray models at 100% power level. The discrete-ordinate radiation 
model discussed in Section 2.3 is used with 12 ordinate directions for these 
simulations. For accurate prediction of radiative heat transfer, the radiative 
properties of combustion gases and walls are required. However, there is a lack of 
knowledge on the radiative properties (absorptivity and emissivity) of MMH and 
NT O. Also, the emissivity and reflectivity of VTE combustor and nozzle walls are 
not known. Due to the lack of this information, the radiative properties for a worst 
case scenario is considered, i.e. the absorptive cross-section of the combustion gas is 
taken to be unity and the wall emissivity is taken to be unity (black). Since the wall 
temperature distribution is not known, the walls are assumed to be maintained at 
1000°K. 
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Figures 5-26 and 5-27 show the temperature distributions predicted using the gas-gas 
and spray models. Comparing these figures with 5-2 and 5-14 shows that the 
temperature distribution is not altered by the inclusion of radiation model. This 
indicates that convective and conductive modes of heat transfer dominate in VTE 
combustion chamber and nozzle. It is found that the source term for the energy 
equation due to radiative heat transfer is an order of magnitude less than those due 
to conduction and convection. This is very clear from Figure 5-28 which compares 
the conductive and radiative fluxes at the combustor and nozzle wall. Figures 5-29 
and 5-30 show the radiative fluxes in the axial and radial directions. It is observed 
that heat is transferred from the combustion chamber to the nozzle region by 
radiation. The predicted ISP for this case is about 314 which is very close the TRW's 
experimental data of 310. The ISP predicted in this case is significantly less than the 
ISPs predicted without the inclusion of radiation model. This is because a constant 
temperature boundary condition is used for all the radiation cases whereas a zero 
gradient boundary condition is used for all the cases with no radiation. The 
conductive and radiative heat loss at the walls leads to a lower ISP. 




Figure 5-26. Temperature Distribution (Gas-Gas Model, 100% Power Level with 
Radiation) 
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Figure 5-27. Temperature Distribution (Spray Model, 100% Power Level with 
Radiation) 
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Figure 5-28. Comparison of Conductive and Radiative Fluxes at Nozzle Wall; o, - 
conductive heat flux; — , radiative heat flux (Gas-Gas Mode, 100% 
Power Level) 
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Figure 5-29. Distribution of Axial Radiation Flux (Gas-Gas Model, 100% Power 
Level) 



Figure 5-30. Distribution of Radial Radiative Flux (Gas-Gas Model, 100% Power 
Level) 
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5.4 Heat Shield Analysis 


The heat shield on VTE has been designed to minimize the radiative heat loss from 
the engine. The amount of radiative heat loss from the chamber and nozzle walls 
depend on wall temperature and their radiative properties. In order to understand 
the effects of heat shield on flow patterns and radiative heat transfer characteristics, 
simulations were performed with and without the inclusion of the radiation model. 

Figure 2-3 showed the 120 X 40 computational grid used for this problem. The 
dimensions for the heat shield were obtained from technical drawings supplied by 
TRW. The grid for VTE is embedded within the grid for shield. The computational 
domain was extended beyond the nozzle exit in order to resolve the expansion of 
nozzle flow into the vacuum environment. Simulations were performed for 100% 
power level with gas-gas combustion model and laminar flow conditions. These 
assumptions are justified since the inclusion of spray dynamics and turbulence is 
not expected to play a role in the heat shield. 

First the flow characteristics inside the heat shield is studied without invoking the 
radiation model. Figure 5-31 shows the Mach number contours for this case. The 
contours inside the combustor and nozzle are very similar to those obtained from 
the calculations without the heat shield. However, the expansion of the gas into the 
heat shield after it leaves the nozzle is observed by the curving of the contour lines. 
The amount of flow entrained into the heat shield is estimated to be negligibly 
small. Figure 5-32 shows the temperature distribution for this case. For one-to-one 
comparison with the radiation results, a constant temperature boundary condition 
is applied for obtaining this figure. A temperature of 1000°K is imposed at the 
chamber/nozzle walls and a temperature of 300°K is imposed at the shield walls. 
Due to this reason, the temperature of the fluid inside the shield region is relatively 

cold, but its density is of the order of 1.0E-05 kg/m 3 . Since the density of the fluid in 
the shield is extremely small, its influence on flow and heat transfer processes are 
expected to be negligible. 
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Figure 5-32. Temperature Distribution (Gas-Gas Model, 100% Power Level, No 
Radiation) 
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For the same problem, the radiative heat transfer characteristics in the shield is 
studied by invoking the radiation model. Since the wall temperature distributions 
are not known, the combustor and nozzle walls are assumed to be maintained at 
1000°K and the shield walls are maintained at 300°K. Figures 5-33 and 5-34 show the 
Mach number and temperature contours for this case. The distributions of Mach 
number and temperature in the chamber and nozzle are very similar to the results 
of the previous case. However, from the Mach number contours, it is observed that 
more fluid is entrained into the shield when radiation effects are included. The 
temperature contours indicate that the plume is relatively cold when compared 
with the predictions without the radiation model. Thus the energy lost due to 
radiation results in a colder plume. Figures 5-35 and 5-36 show contours of radiative 
fluxes in the axial and radial directions. The heat radiated by nozzle wall into the 
shield region can be clearly seen. It can also be seen that heat is being radiated away 
from the combustor region into the nozzle. Since the radiative heat flux at the exit 
is positive, some amount of heat from the nozzle region is lost by way of radiation. 

Figures 5-37 and 5-38 compare the conductive and radiative heat fluxes at the nozzle 
and outer shield walls. At the nozzle wall, the convective heat fluxes are much 
higher than the radiative heat flux. On an average, the radiative heat transfer is 
found to be less than 8% of the conductive heat transfer. On the other hand, the 
radiative heat flux at the outer shield wall is much higher than the conductive heat 
flux, but small when compared with the heat flux at the nozzle wall. Thus small 
amount of heat is lost from the shield to outer space due to radiation. 
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Figure 5-33. Mach Number Distribution in VTE (Gas-Gas Model, 100% Power Level, 
with Radiation) 
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Figure 5-34. Temperature Distribution (Gas-Gas Model, 100% Power Level, with 
Radiation) 
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Figure 5-35. Distribution of Axial Radiative Flux (Gas-Gas Model, 100% Power 
Level) 
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Figure 5-36. Distribution of Radial Radiative Heat Flux (Gas-Gas Model, 100% 
Power Level) 
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Figure 5-38. Comparison of Conductive and Radiative Fluxes at Shield Wall; 

o, - conductive flux; — , radiative flux (Gas-Gas Model, 100% Power 

Level) 


69 


4075/3 


6. CONCLUSIONS AND RECOMMENDATIONS 
6.1 Accomplishments 

a. A two-liquid, two-phase Eulerian-Lagrangian spray model has been 
incorporated into the REFLEQS code. 

b. The spray tracking methodology has been extended to account for complex, 
non-orthogonal BFC grids and interaction of the spray with the domain 
boundaries. 

c. A complete coupling between the gas phase and the liquid droplets has been 
implemented through source/sink terms for the mass, momentum, and 
energy equations of both phases. 

d. A chemical equilibrium model (based on the element potential method) has 
been incorporated into the code. This model has been validated against the 
predictions of the TDK code. 

e. A linear stability analysis model has been included to estimate drop sizes 
resulting from the impingement of two axisymmetric liquid sheets. 

f. A stochastic droplet dispersion model has been implemented into the code to 
model the droplet-turbulence interaction. 

g. A hot gas radiation model (based on the discrete ordinate approach) has been 
incorporated into the code. The model has been validated for simple 
problems having exact solutions. 

The above capabilities in conjunction with the existing capabilities of REFLEQS 
have facilitated the simulation of complex physical phenomenon in liquid rocket 
engines such as OMV/VTE combustors. However, certain limitations and 
shortcomings do exist that need to be addressed in the future. 


70 


4075/3 


6.2 Summary of Findings from the Numerical Analysis 


A number of simulations were obtained to analyze the performance of VTE at 
different power levels. Predictions using the gas-gas and the liquid-gas spray models 
have been compared. In most of the cases (both gas-gas and spray models), the 
predicted ISP is approximately 320 which is slightly higher than the measured value 
of 310. The ISP is overpredicted because the heat loss from the engine walls is not 
included in our study when adiabatic boundary conditions are imposed at the walls. 
In general, the gas-gas model predicts a higher ISP than that predicted by the spray 
model, because the energy required to atomize and evaporate the droplets are not 
present in the gas-gas model. For both the models, the predicted ISP slightly 
decreases at off-design power levels indicating a decrease in efficiency. 

The predictions of spray model suggest that smaller droplets evaporate faster, 
leading to better mixing in the combustion chamber. This results in better 
performance of the engine. The larger droplets, on the other hand, are shown to 
impinge on the splash plate providing a cooling effect for the wall. It is shown that 
optimum performance can be achieved when the droplet sizes are in the 25-50 
microns range. The predictions of gas-gas model show two recirculation zones: one 
in the swept bach region of the combustor and the other in the center of the 
chamber. These recirculation zones are absent in the predictions of the spray model. 
The recirculating zones in the predictions of gas-gas model are mainly induced by 
the relatively high inlet velocity of the gases. The predictions also show that the 
hottest region is at the center of the combustion chamber. As a result, the head end 
of the pintle injector is expected to be very hot. This has been observed in VTE hot 
fire tests. Later designs of the VTE engine has an ablative material covering the 
head end. 

Predictions using the radiation model suggest that radiation effects on flow and heat 
transfer inside the VTE combustion chamber and nozzle are insignificant. The 
radiative heat flux at the walls are an order of magnitude less than the conductive 
heat flux. The results show that heat is transferred from the combustion chamber 
zone to the nozzle region by radiation. Simulations performed with the heat shield 
show that negligible amount of fluid is entrained into the heat shield region. 
However, the heat shield is shown to be effective in protecting the OMV structure 
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around the engine from the radiated heat. Due to the limitted knowledge of the 
radiative properties of the OMV combustion gases, very high values of gas 
absorptivity and wall emissivity are used. For more realistic values of radiative 
properties, the contribution due to radiative heat transfer is expected to be smaller. 

6.3 Su ggestions for Future Work 

a. An advanced model for atomization due to impingement of unlike liquids in 
sheet (2-D) and round jets (3-D) forms. In the present analysis, a simplified 
surface wave model was used to estimate the mean droplet size. This surface 
wave model considers the instability of waves on like impinging jets in a 
stagnant environment. On the other hand, in the VTE engine, unlike jets 
impinge in a reacting gas environment. Hence the present simplified 
atomization model overpredicts the droplet size. A rigorous model including 
the effects of ambient reacting gas is needed in order to accurately predict 
the droplet size. 

b. The development of a coupled atomization model . In the actual conditions 
inside the VTE engine, the droplet sizes are expected to be a strong fuction of 
local gas conditions. The relative velocity between the liquid sheet and the 
gas will vary along the sheet radius resulting in droplet size variation. The 
simplified atomization model, on the other hand, only predicts a single mean 
size based on the most unstable wavelength. Thus, a technique such as the 
Jet-Embedding technique 23 ' needs to be developed for atomization in 
impinging jets. 

c. Development of advanced models to simulate mixing and combustion in the 
liquid phase. When unlike jets impinge on each other, the resulting droplets 
are inhomogeneous. In order to predict the composition of droplets, the 
mixing occuring in the liquid phase needs to be modeled. This requires a 
solution of the liquid phase equations near the impingement point. Since the 
liquid phase equations near the impingement point are strongly elliptic and 
since free surfaces are involved, obtaining a solution for this problem would 
represent a major effort. 
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d. 


Identification and implementation of reduced finite rate reaction 
mechanisms for complex fuel /oxidizer combinations. In the calculations of 
the present study, an equilibrium chemistry model was used to simulate 
combustion. For steady state calculations, this model was shown to provide 
satisfactory results. However, for studying transient phenomena such as 
combustion instabilities, a finite rate reaction models are necessary. For the 
reaction between MMH and NTO, several reaction steps with 200 species are 
involved. A detailed modeling of this reaction would be prohibitively 
expensive. Hence reduced finite rate reaction models need to be developed 
for MMH-NTO and other commonly used propellants such as H2-02. 

e. Implementation of dense spray models to simulate the spray in the region 
close to the injector. The present spray model assumes that the volume 
occupied by the droplets are small when compared to the cell volume. This 
assumption does not hold good near the injector where the spray is dense. 
The dense spray effects are thus need to be included in the calculation 
procedure. 

f. Development and incorporation of models to account for supercritical 
behavior of fluids. Engines such as SSME operate under supercritical 
pressure conditions. The effects of supercritical conditions on evaporation 
and heat transfer are shown to be important due to the abrupt changes in 
fluid properties. Several models are available for supercritical vaporization 
and heat transfer. These models should be incorporated in the present code. 

g. Extend the droplet dispersion model to non-isotropic turbulence. The 
turbulent droplet dispersion model currently implemented in the code 
assumes that turbulence is locally isotropic. This assumption is strictly not 
valid in the VTE combustion chamber. For accurate prediction of droplet 
dispersion, the effect of non-isotropy must be considered. 
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Appendix A 


Sample Input Files 


auto2d 

* Input data for the OMV Variable Thrust Engine 

******* ‘********Q as _(3. as ^ 0< j e ]**=4c*>K***^** ***************** 

title OMV-VTE (100x30,modvil =3, turbulent, MMH/N204) 

* Parameters * 

* ox stands for oxidizer, fu stands for fuel, 
integer L=100 

integer M=30 

real uox= 110.52 vox=-51.53 

real ufu=0.0 vfu= 144.3 

real Tiniet=300.0 Pinlet=6.89250E+05 

* Geometry & Grid * 

1 L ; m M 

polar T ; bfc T ; ortho g F 

* 

grid from omvGRID 

* Solution control * 

niter 1000 ; info 0 

algorithm simplec ; scheme upwind all 

* 

solve nvu nw nvpp nvkl nvdl nvh 
* 

absor 1.0; scatr 0.0; emisw 1.0 
* 

niterp 5 ; iblox nvpp 0 ; ibloy nvpp 0 
* 

nsweep nvu 10 ; nsweep nw 10 ; nsweep nvpp 30 ; nsweep nvkl 10 

nsweep nvdl 10 ; nsweep nvh 10 ; nsweep nvfl 10 ; nsweep nvf2 10 
* 

dtfals nvu 1.5 ; dtfals nw 1.5 ; dtfals nvkl 1.5 ; dtfals nvdl 1.5 

dtfals nvh 1.5 ; dtfals nvfl 1.5 ; dtfals nvf2 1.5 
* 

relax nat 0.2 ; relax narho 0.2 ; relax nap 0.2 ;relax navis 0.2 
fimax navis 0.1 

* Output Control * 

Iprint nvu 0 ; Iprint nw 0 ; Iprint nap 0 
cuser T ; pltdmp t ; Idump 50 

ildbg 50 ; ildbg 56 ; jldbg 25 ; jldbg 30 

* Physical properties and models * 

press 0.E+05 ; moden 3 ; comprs T ; mixing 2 ; modvds 0 
modvis 2 ; modwal 3 ; modvil 3 ; ICP 3 ; iequil 1 ; modken 1 

* * 

ncmps 2 

compos 1 n2o4 1.0 h 0.0 h2 0.0 h2o 0.0 no 0.0 co 0.0 ch4 0.0 
compos 2 ch6n2 1.0 n2 0.0 o 0.0 oh 0.0 o2 0.0 co2 0.0 
* Initial Conditions * 
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uin 100. ; vin 50. ; presin Pinlet ; temin 2000. ; compin 7 
♦restart from MODEL.AUR 

* Boundary region setup * 

nblock 1 

poros vol 0.0 1 7 1 5 

* * 

nreg 9 

region 1 west fixm 116 9 

region 2 west wall 1 1 10 30 

region 3 west wall 8 8 15 

region 4 east extrap L L 1 M 

region 5 south symm 8 L 1 1 

region 6 south wall 1 1 6 6 

region 7 south fixm 2 4 6 6 

region 8 south wall 5 7 6 6 

region 9 north wall 1 L M M 

♦ Boundary Condtion Specifications * 

uvela uox ; wela vox ; presa Pinlet ; tema Tinlet ; compa 1 

uvelb ufu ; welb vfu ; presb Pinlet ; temb Tinlet ; compb 2 

be vu = va 

be set 1 vu 

be vu = vb 

be set 7 vu 

be vu = 0.0 

be vu nvh = zgbc 

be set 2 vu 

be set 3 vu 

be set 4 vu 

be set 5 vu 

be set 6 vu 

be set 8 vu 

be set 9 vu 

end 
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auto2d 

* Input data for the OMV Variable Thrust Engine 
****** * ******2prQy Model********** ****** ******** 
title OMV-VTE ( 100x30, modvil =3, turbulent, MMH/N204) 

* Parameters * 

* ox stands for oxidizer, fu stands for fuel, 
integer L=100 

integer M=30 

real Tinlet=300.0 Pinlet=6.89250E+05 

* Geometry & Grid * 

I L ; m M 

polar T ; bfc T ; orthog F 

* 

grid from omvGRID unform 

* Solution control * 

niter 1000 ; info 0 

algorithm simplec ; scheme upwind all 

* 

solve nvu nw nvpp nvkl nvdl nvh 
* 

niterp 5 ; iblox nvpp 0 ; ibloy nvpp 0 
* 

nsweep nvu 10 ; nsweep nw 10 ; nsweep nvpp 30 ; nsweep nvkl 10 

nsweep nvdl 10 ; nsweep nvh 10 ; nsweep nvfl 10 ; nsweep nvf2 10 
* 

dtfals nvu 1.0 ; dtfals nw 1.0 ; dtfals nvkl 1.0 ; dtfals nvdl 1.0 

dtfals nvh 1.0 ; dtfals nvfl 1.0 ; dtfals nvf2 1.0 
* 

relax nat 0.3 ; relax narho 0.3 ; relax nap 0.3 

* Output Control * 

cuser T ; pltdmp t ; Idump 50 

Iprint nvu 0 ; Iprint nw 0 ; Iprint nap 0 

* Physical properties and models * 

press 0.E+05 ; moden 3 ; comprs T ; mixing 2 ; modvds 0 
modvis 2 ; modwal 3 ; modvil 3 ; ICP 3 ; iequil 1; modken 1 

* * 

ncmps 2 

compos 1 n2o4 1.0 h 0.0 h2 0.0 h2o 0.0 no 0.0 co 0.0 ch4 0.0 
compos 2 ch6n2 1.0 n2 0.0 o 0.0 oh 0.0 o2 0.0 co2 0.0 

* Initial Conditions * 

uin 100. ; vin 50. ; presin Pinlet ; temin 3000. ; compin 7 
*restart from MODEL.AUR 

* Spray data * 

spray t 

tboill 300. ; cpll 1580.0 ; hill 414000.0 ; rholl 1433.0 
psatll 1.2E5 


A-4 


4075/3 


tboil2 360. ; cpJ2 3000.0 ; hJ12 875000.0 ; rhol2 870.0 

psatl2 6.7E3 

iuser 52 1 ; iuser 53 5 

nslot 2 

* slot i * 

inject slot 1 20 1 0.1460 CELL 4 7 

inject size 1 lognormal 50.0e-6 25.0e-6 100.0e-6 0.2 

inject be 1 25.0 27.0 0.0 300.0 

* slot 2 * 

inject slot 2 20 2 0.0890 CELL 4 7 

inject size 2 lognormal 50.0e-6 25.0e-6 100.0e-6 0.2 

inject be 2 25.0 23.0 0.0 300.0 

* Boundary region setup * 

nblock 1 

poros vol 0.0 1 7 1 5 


nr eg 10 
region 1 
region 2 
region 3 
region 4 
region 5 
region 6 
region 7 
region 8 
region 9 


west 

west 

west 

west 

east 

south 

south 

south 

south 


wall 

fixm 

wall 

wall 

extrap 

svmm 

wall 

fixm 
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region 10 north wall 

♦ Boundary Condtion Specifications * 

uvela 0.0 ; wela 0.0 ; presa Pinlet ; tema Tinlet ; compa 1 

uvelb 0.0 ; welb 0.0 ; presb Pinlet ; temb Tinlet ; compb 2 

be vu = va 

be set 2 vu 

be vu = vb 

be set 8 vu 

be vu = 0.0 

be set 5 vu 

be vu = 0.0 


be vu nvh = zgbc 
be set 1 vu 
be set 3 vu 
be set 4 vu 
be set 6 vu 
be set 7 vu 
be set 9 vu 
be set 10 vu 
end 
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